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Abstract 

We examine the possibility that a metastable quantum state could ex- 
periment a phenomenon similar to thermal activation but at zero temper- 
ature. In order to do that we study the real-time dynamics of the reduced 
Wigner function in a simple open quantum system: an anharmonic oscilla- 
tor with a cubic potential linearly interacting with an environment of har- 
monic oscillators. Our results suggest that this activation-like phenomenon 
exists indeed as a consequence of the fluctuations induced by the environ- 
ment and that its associated decay rate is comparable to the tunneling rate 
as computed by the instanton method, at least for the particular potential 
of the system and the distribution of frequencies for the environment con- 
sidered in this paper. However, we are not able to properly deal with the 
term which leads to tunneling in closed quantum systems, and a definite 
conclusion cannot be reached until tunneling and activation-like effects are 
considered simultaneously. 

1 Introduction 



The study of the decay rate of a state trapped in a metastable state by a poten- 
tial barrier has a long and distinguished history both in statistical physics and 
in quantum mechanics. In statistical mechanics one is usually worried about the 
thermal activation effect, by which a particle escapes over a potential barrier 
due to the fluctuations induced by a thermal bath. As a paradigm we have the 
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classical work by Kramers (1940), who considered a classical Brownian parti- 
cle trapped in a metastable minimum and computed the escape probability by 
analyzing the dynamics of its probability distribution function, both in the un- 
derdamped and overdamped cases. In quantum mechanics one is interested in 
the tunneling effect, by which a particle escapes the local minimum by traversing 
the potential barrier through a classically forbidden region. A technique which 
has been of great success to compute the tunneling rate in quantum mechanics 
and quantum field theory is the instanton method (Coleman, 1977; Callan and 
Coleman, 1977; Coleman, 1985), where the decay probability can be computed in 
terms of the classical trajectories in imaginary time. In both thermal activation 
and tunneling the decay rate r, meaning the decay probability per unit time, 
follows an approximate exponential law, r ~ Aexp (—B). In thermal activation, 
-B = Vs/{kT), where k is Boltzmann's constant, T is the absolute temperature 
and Vg is the height of the free energy measured from the metastable minimum; 
in tunneling B = Se/^, where Se is the action for a suitable trajectory which 
goes under the barrier in imaginary time. 

In recent years mesoscopic physics has become a center of attention. Ex- 
perimental advances are pushing the boundaries between classical and quantum 
systems and lead, in particular, to the possibility of observing quantum tun- 
neling for systems that can be described by macroscopic variables. These are 
essentially open quantum systems, which are characterized by a distinguished 
subsystem within a larger closed quantum system, described by some degrees 
of freedom which are subject to physical experimentation, and the rest of the 
system, described by generally unobservable degrees of freedom which act as an 
environment or bath for the distinguished subsystem. The environment induces 
both dissipation and noise to the distinguished system, which is usually referred 
to as the "system" for short. Many of these open quantum systems can be equiv- 
alently represented by a particle subject to an arbitrary external potential and 
coupled to an environment consisting of an infinite set of independent harmonic 
oscillators. A number of known physical systems can be modelled by adjusting 
the coupling of the system and environment variables and choosing appropriate 
potentials. 

Caldeira and Leggett (1981, 1983b), in two influential papers, considered the 
effect of an environment on quantum tunneling. They were able to generalize 
the instanton method to a simple open quantum system. In particular, they 
considered an anharmonic quadratic plus cubic potential bilinearly coupled an 
environment at zero temperature consisting of an infinite set of harmonic oscilla- 
tors, with frequencies distributed according to the so-called ohmic distribution. 
They argued that this system is a very good model for the flux trapped in 
a superconducting quantum interference device (SQUID), a single Josephson 
junction biased by a fixed external current, and others. Assuming that the envi- 
ronment degrees of freedom are only weakly perturbed by the interaction with 
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the system, they concluded that dissipation always tends to suppress tunneling. 

Fujikawa et al. (1992a,b) reanalyzed the same problem as Caldeira and 
Leggett employing canonical perturbation theory and the second quantization 
formalism for the system instead of the instanton method, with a quadratic plus 
quartic potential in order to have a well-defined ground state. By considering 
the two lowest energy eigenstates of the system they were able to reproduce 
Caldeira and Leggett's results. They studied the effect of the next two excited 
states and concluded that, once these states are taken into account, dissipation 
can enhance tunneling for some distributions of frequencies of the environment, 
but not for the ohmic distribution considered by Caldeira and Leggett. 

In some recent papers a different real-time approach to compute the vacuum 
decay rate in quantum field theory was introduced (Calzetta et al, 2001, 2002) 
as a first step to consider situations far from the equilibrium. The analysis is 
based on the time evolution equation, the so-called master equation, for the 
reduced Wigner function describing the open quantum system. The quantum 
field theory problem was reduced to an open quantum system described by a 
single degree of freedom associate to the modes of the field which are nearly 
homogeneous within a region whose size corresponds to a nucleating bubble of 
true vacuum, and coupled to an infinite set of harmonic oscillators corresponding 
to the inhomogeneous modes of the field. The coupling was linear in the system 
variables but quadratic in the environment ones. The master equation contained 
dissipation and noise terms, which describe the influence of the environment on 
the system, as well as derivative terms in the momentum coordinates (third 
order for a cubic potential), which are responsible for quantum tunneling in 
closed quantum systems. 

Unfortunately it is not possible to compute the total decay rate in a closed 
analytic or semianalytic form because it is very difficult to deal with the third 
order derivative terms. One should resort to methods such as those based on 
matrix continued fractions in order to compute the decay rate master equations 
with third-order derivative terms (Risken and Vogel, 1988). However in the work 
by Calzetta et al. (2001, 2002), as well as in the present paper, one was primarily 
interested in the contribution from the environment backreaction on the vacuum 
decay rate and, hence, the rather drastic approximation of neglecting the third 
derivative term responsible for tunneling was made. The decay rate obtained 
was entirely due to the terms which at high temperature are responsible for 
thermal activation. In principle this approximation should be correct provided 
that the timescales associated to the activation-like effect and the tunneling 
effect are very different. The remarkable result was that the activation-like 
effect produced by the backreaction of the inhomogeneous modes was, in fact, 
larger that the tunneling effect obtained with the instanton method. 

The fact that activation may be important as a backreaction effect can be 
understood by noting that the characteristic timescale for the decay process is 
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much larger than the dynamical and relaxation timescales; see Calzetta and 
Verdaguer (1999) for a detailed analysis. In fact, although the effects of dissi- 
pation and noise are very small on the characteristic dynamical timescale, they 
can have a cumulative effect which becomes important in the long run. The real 
time approach based on the master equation seems a suitable technique to deal 
with those backreaction effects, but is difficult to implement when addressing 
the tunneling effect. On the other hand, the instanton method is very well suited 
for studying the tunneling effect, but, at least in some cases, seems to underplay 
the backreaction of the environment on the system. 

The aim of this paper is to study whether the results obtained by Calzetta 
et al. (2001, 2002) are a particular feature of the particular system that they 
considered or, on the contrary, a general feature which can be extended to a 
wider class of systems. Therefore, using the real-time techniques mentioned 
above, we reanalyze the simpler model studied by Caldeira and Leggett (1981, 
1983b), in order to check whether a similar activation- like effect is also found in 
that case. 

There are several differences between the model considered in this paper 
and the field theory case analyzed by Calzetta et al. (2001, 2002). First, we 
consider here a bilinear coupling between the system and environment, whereas 
the coupling considered in the field theory case was quadratic in the environment 
degrees of freedom. Second, in the field theory case both the spectral distribution 
of the environment frequencies and the value of the coupling constant were a 
priori determined by the particular system-environment separation considered 
there, whereas we shall freely choose the spectral distribution and the coupling 
parameter. Third, we will employ several techniques (harmonic approximation, 
Kramers method and lowest eigenvalue expansion) which could not be applied 
to the field theory problem because of the particular features of that model. 

The plan of the paper is the following. In section 2 we set up the model. 
In section 3, we present the master equation for the reduced Wigner function. 
Then, following closely Calzetta et al. (2002), we neglect the third order deriva- 
tive term and concentrate on the weak dissipation case by studying the averaged 
dynamics over an oscillation period. In section 4 we obtain an analytical expres- 
sion for the decay rate by assuming a harmonic approximation for the classical 
trajectories. Two alternative approaches are employed. The first one is based 
on a perturbative expansion for the lowest eigenvalue, whereas the second one 
is based on Kramers's classical work. Finally, in section 5 we discuss the impli- 
cations of our result. 

2 The open quantum system model 

Let us consider a particle of mass M, the "system", subject to an arbitrary 
potential V{x) and coupled to a bath of independent harmonic oscillators of 
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mass m, the "environment". Let us assume that the system and the environment 
are Unearly coupled. The action for the whole set of degrees of freedom is given 

by 

{qj}] = Ss[x] + Sc[{qj}] + Sint[x, {qj}], (la) 

where the terms on the right-hand side, which correspond to the action of the 
system, the environment and the interaction term respectively, are given by 

Ss[x] = Jdt Qmx^ - Vix)^ , (lb) 

Seiiqj}] = E/d^ [Imq] - ^mu;|g|] , (Ic) 
j 

Sint[x,{qj}] = ^Cj / dtx{t)qj{t), 

with Cj being the system-environment coupling parameters and Uj the environ- 
ment oscillator frequencies. The system potential V{x) includes a quadratic 
part, corresponding to an oscillator of frequency fio, and an anharmonic part 

y{nl)^ 

V{x) = hllx^ + V^^'\x). (2) 

At this point, the potential V^^^\x) is arbitrary, but later on we will take a 
cubic potential V^^^\x) = — (A/6)x^. It will be convenient for us to rewrite the 
interaction term as 

SiAx,{qj}]= [ ^^I{io) I dtx{t)qit;u;), (3) 
Jq ■kc{uj) j 

where c{lo) and q{t;Lo) are functions such that c{LOj) = Cj and q{t;LOj) = qj{t) 
and 

j ■' 

is the spectral density of the environment. 

When the system and the environment are initially uncorrelated, i.e., when 
the initial density matrix factorizes, the evolution for the reduced density matrix 
can be written as 

(^X ^ X J ^ ) — ^ diC| J {^X ^ X ji/'f X\ ^ X'y ^ Pf ("^i 5 X'y j j 

where the propagator J is found to be, in a path integral representation, 

xit) = Xf x'{t) = x'f 

J{xf,x'i,t;xi,x{,tO = J Vx J Pa;'e'(^M-^[^']+^i''[^'^'])/'^ ^6) 

x{ti)=Xi x'{ti)=x[ 
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where -Sipfx, x'] is the influence action, related to the influence functional Fjp in- 
troduced by Feynman and Vernon (1963) through Fip[x, x'] = exp {iSiF[x, x']/h). 
For a Gaussian initial density matrix for the environment, the influence action 
can be expressed as (Feynman and Vernon, 1963; Feynman and Hibbs, 1965; 
Caldeira and Leggett, 1983a): 



Sif[x,x']=~2 f ds j ds'A(s)L>(s,s')S(s') 
+ '-j\s j\s'/^{s)N{s,s')/^{s'), 



where S = {x + x')/2 and A = x' — x. The kernels D{t, t') and N{t, t') are called 
the dissipation and noise kernels, respectively. If initially we are at thermal 
equilibrium at a temperature T these kernels are given by: 

D{t,t')= I{u;)smu;{t-t'), (8a) 

Jo ^ 

N{t,t') = ^/(o;) coth j cosa;(t - t'). (8b) 

The influence action can be divergent and a rcnormalization procedure may 
be required, as can be seen by reexpressing the influence action as 



Sw[x,x'\ = j ds [ ds' A{s)H{s,s')J:{s') 
Jti Jti 

ds / ds' A{s)N{s,s')A{s'), 



(9) 



where, at least formally, H{t,t') = -29(1 - t')D{t,t'), being e{t - f) the step 
function. The kernel H{t, t') is a product of two distributions, which in general is 
not well defined and may contain divergences. Nevertheless, it is always possible 
to introduce suitable counterterms in the bare frequency of the system JIq in 
order to compensate the divergent terms coming from i?(t, t'). See Caldeira and 
Leggett (1983a); Roura and Verdaguer (1999); Calzetta et al. (2003) for more 
details. However in the particular problem in which we are interested this issue 
will turn to be unimportant since the divergent parts of the kernel H{t, t') will 
cancel in the final results. Thus, we can use the bare kernel H{t,t') (with some 
implicit regularization) instead of its renormalized expression. 

Following Caldeira and Leggett (1981, 1983b) we shall consider the case of 
zero temperature and ohmic environment, in which we have a continuum of 
harmonic oscillators in the environment distributed according to: 

I{lo) = r]uj. (10) 
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With this spectral density the expectation value of x{t) obeys the equation of 
motion of a classical damped oscillator with a friction coefficient given by the 
proportionality constant rj. In this case dissipation and noise kernels are found 
to be: 



D{t,t') = r)S'{t-t'), (11a) 

TT (t - 



^(«'*') = zPf7^^> (lib) 



where Pf indicates the Hadamard finite part prescription (Schwartz, 1957). 
Later on we will need the expressions of the Fourier transforms of the noise 
and dissipation kernels, 

D{t,t') = I ^e— (t-0^(a;), (12a) 
" / ^ e-^'^(*-*') N{lo), (12b) 
which in the case of zero temperature and ohmic environment are given by 

D{(jj) = irioo, N{uj) = r?|w|. (13) 



V 




Figure 1: Plot of the potential V{x) under which the particle is confined. The 
maximum of the potential barrier is at Xg and corresponds to an escape energy 
Es, which we consider much larger than the zero point energy of the harmonic 
oscillator, Klo/2. 

We will concentrate on the potential V{x) = {l/2)M^'^x'^ — (A/6)a;^, which 
exhibits a metastablc minimum at ,t = and an unstable maximum at x = = 
2Mf2^/A, which corresponds to an energy = V(xs) = 2M^ri^/(3A^) (see Fig. 
1). The system state will be peaked located around the metastable minimum at 
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X = 0, and can escape through the potential barrier. We will consider that once 
the particle exits the potential well region, it never reenters. Since this potential 
is not bounded from below, it should be understood as an approximation to a 
more realistic situation in which there exists an absolute minimum, located at a 
much lower energy, so that the return probability is negligible. We will restrict 
to the situation in which the energy barrier is much larger than MIq. 

3 Phase-space dynamics 

3.1 Evolution of the reduced Wigner function 

The reduced Wigner function is a phase space distribution defined from the 
reduced density matrix pr by the following integral transform: 



The Wigner function is a quantum mechanical analogue of a phase space proba- 
bility distribution (Wigner, 1932; Hillery et al, 1984). The partial distribution 
dpW,:{x,p,t) gives the probability density of finding the system at the po- 
sition x; in the same way, f^^dxWr{x,p,t) gives the probability density of 
finding the system with momentum p. However, the uncertainty principle pre- 
vents us from determining at the same time the position and the momentum of a 
particle, so that the Wigner function cannot be interpreted as a true probability 
density in phase space. In fact, the Wigner function it is not necessarily positive 
defined everywhere and in general it may acquire negative values. 

Up to terms of order cj, he?', and with Cj being the environment coupling 
constants, the reduced Wigner function Wj- = Wj.{x,p, t) obeys the following evo- 
lution equation (Calzetta et al, 2002; Roura, 2001), which we shall call master 
equation: 




(14) 



at 



24 ap^ ' 



(15) 



where {•, - jpH are the Poisson brackets. 



{f,9} 




Hg is the reduced system Hamiltonian, 




(16) 
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and T> and M are given by 

V{t) = -2 j^^ dt' D{t,t')X{t';x,p), (17a) 

M{t)= [ dt' N{t,t')X{t';x,p), (17b) 

where X{t';x,p) is a solution of the classical equations of motion associated to 
the hamiltonian Hg, 

dX P , , 

-dF=M^ ^''^^ 

^ = -MnlX + A—, (18b) 

with final conditions X{t) = x and P{t) = p. 

For the case of ohmic environment the friction coefficient ij is proportional 
to Cj, and the above approximation for the evolution equation will be valid and 
consistent when the condition S ^ h > S^y /Qq is fulfilled, where 7 = ri/{2M) 
is the characteristic dissipation frequency and S is the typical action for the 
process considered. If we are interested in studying the jump over a potential 
barrier of height Eg this condition can be restated as Eg ^ hO,Q > £sj/0,q, which 
shows that the master equation will be valid for small dissipation and large 
energy barriers. 

At this point it is worth making a comment on the notation. Throughout 
this paper we will use lowercase letters {x,p, 0,j.. .) to indicate phase-space vari- 
ables, which are the arguments of phase-space distributions such as the Wigner 
function, whereas the corresponding uppercase letters {X, P, 6, J . . .) will indi- 
cate time-functions which give the phase-space position of a particle in a given 
time. 

In order to introduce the Fourier transform of the coefficients V and J\f later 
on, it will be convenient for us to replace the integration limit in Eqs. (17a) 
and (17b) by —00, although the initial conditions are set up at ti = 0. For the 
coefficient T> this does not introduce any error, since the dissipation kernel only 
has support in t = t'. For the coefficient M, this introduces a small error, which 
can be estimated by performing the integral in Eq. (17a), choosing a periodic 
function for X{t') of frequency Oq; this will be enough in the context of the 
adiabatic approximation, which we will introduce later on. The result of the 
calculation shows that for times t which verify t 3> ^ the contribution of 
the integral from —00 to is comparatively small. Since at the end we shall be 
interested in studying time scales of the order of the decay time, which are much 
larger than the characteristic dissipation time 7"^, which will in turn be much 
larger than the time of oscillation ^, the approximation can be considered 
safe. 
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If the system were isolated, Eq. (15) would reduce to 

where W is the Wigner function of the closed system. This equation is exactly 
equivalent to von Neumann's equation for the density matrix of a one dimen- 
sional quantum mechanical system with the potential V{x). If the last term in 
this equation were not present, the evolution of the Wigner function would be 
entirely equivalent to that of a classical ensemble in the phase space. Hence, 
the term with the third derivatives must be the responsible for tunneling. In 
principle, one could compute the tunneling amplitude from Eq. (19), but in 
practice tunneling is more easily calculated with the WKB approximation to 
the Schrodinger equation (see, e.g., Galindo and Pascual, 1991; Landau and 
Lifshitz, 1977) or the instanton method (Coleman, 1977; Callan and Coleman, 
1977; Coleman, 1985). As we have discussed, in this contribution we are not 
going to deal with tunneling, but rather to compute the effect due to activation. 
With this aim, we shall neglect the last term of Eq. (15), which is the responsi- 
ble for tunneling. Hence we will use the following equation for the distribution 
function Wr'. 



= {H„ WjpB + |- {VW, + h{M, H^rjpB) . (20) 



Formally Eq. (20) can be thought as a Fokker-Planck equation, describing 

the dynamics of an ensemble of points in the phase space. The dynamics of 
this ensemble of points can be equally characterized by means of the following 
Langevin equation: 

X. = (21a) 

/oo 
dt'Hit,t')X,{t') + ^ 

-~ (21b) 

where ^ is a Gaussian noise of zero mean and correlation function 

{amt'))^ = hN{t,t'), (21c) 

Xg = Xs{t, ^] and Pg = Ps{t, ^] are the stochastic functions corresponding to the 
phase space variables x and p respectively, and the second equality in Eq. (21b) 
is just valid in the case of ohmic environment. This Langevin equation docs 
not describe actual trajectories of the system (meaning a continuous sequence 
of projectors for the position and the momentum of the system at each instant 
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of time, which would violate Heisenberg's uncertainty principle) but must be 
regarded as a formal computational tool, in the same way as the Wigner func- 
tion does not correspond to a true probability density. In fact, Langevin-like 
equations appear naturally in the context of open quantum systems when trying 
to derive the dynamics of the reduced Wigner function from the path integrals 
in Eq. (5) (Calzetta et al, 2003). 



3.2 Action-angle variables 

In order to obtain explicit expressions for the coefficients V and TV, wc need to 
solve the set of Eqs. (18), which describe the motion of a classical particle in the 
potential V{x). Since we are interested in the motion inside the potential well, 
which is periodic, it is possible to introduce action-angle variables 9 and j (see, 
e.g., Goldstein, 1980). The action variable is defined by 



IfpdX, (22) 



•' 27r 

whereas the angle variable 6 changes from zero to 27r and is canonically conjugate 
to j. Recall that the reduced system Hamiltonian can be entirely written in 
terms of the action variable, Hs = Hs{j). 

We shall consider 6 and j phase-space variables like x and p, and we will 
analyze the Fokker-Planck equation in terms of these new variables. However, 
the solution of Eqs. (18) can be also described by giving the trajectory of the 
particle in the 6 j space. In a completely analogoTis way to X and P, we will 
consider the functions Q{t';9,j) and J{t';9,j), which give the angular position 
and action of a particle satisfying the set of Eqs. (18), with final conditions 
0{t) = 9 and J{t) = j. The functions Q and J satisfy the equations of motion 
= io{J) and J = 0, where ui{j) = dHs{j)/dj is the frequency of oscillation. 
With the aforementioned boundary conditions, these equations of motion can 
be immediately solved to give Q{t'; 9,j) = 9 + fl{j){t' — t) and J{t'; 9,j) = j. 

Since 9 is an angle, the transformation equation x = x{9,j) is periodic in 9, 
x{9,j) = x{9 + 27r,j) and thus it can be decomposed in terms of a Fourier series 
with respect to 9, 

oo 

x{9,j)= J2 e'"'^n(i), (23) 

n=— oo 

where X-n{j) = x'^ij) since x is real. The trajectory of the particle can be also 
decomposed in terms of the Fourier series associated to the angular coordinate: 

oo 

X{t';9,j)= e^^^(''MxniJ{t';9,j)) 

(24) 

^ J2 e^"[^+^o(*'-*)la;„(j). 

n=— oo 
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Then we can write the functions V{t) and Af{t) appearing in Eq. (20) as 



OO 



dt' D{t,t')X{t') = J2 



(25a) 



— oo 



n=— oo 




(25b) 



where Dn{j) and Nn{j) are given by 




dio 2iD{uj) 



(26a) 



n 



27r o; + n^l{j) — ie' 
duj -iN{uj) 



n 



2iT uj + nQ{j) — ie 



(26b) 



In order to derive these last expressions we have made use of the fohowing 
equahty: exTp{isu) = i/{s + ie), where e is an arbitrarily small positive real 
number. 

3.3 Weak dissipation limit: averaging over angles 

As discussed by Kramers (1940), in the case of small dissipation, i.e., 7 ^ Qq, 
the phase space dynamics will mostly correspond to a gradual change of the 
distribution of the ensemble over the different energy values. The change of 
the Wigner function over an oscillation period will be small, so that we may 
suppose that the reduced Wigner function only depends on the action variable 
j (or the energy E), and does not depend on the angular variable 6, Wr{6,j) = 
F{j). Thus, we can obtain a simpler equation by averaging all the terms of 
the Fokker-Planck equation over the variable 6. Furthermore, in this case the 
averaged Wigner function F is a partial distribution, and hence it admits a true 
probabilistic interpretation, as opposed to the non-averaged Wigner function. 
Notice that in this case {Hg, F}pb = and that, for any phase space function 



Using the last expression we can average Eq. (20) to obtain the following aver- 
aged Fokker-Planck equation: 





(27) 
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where we have introduced ^^{j) and which are defined as fohows: 

oo 

i J2 K{j)\^nDn{j), (28a) 

n=—oo 

oo 

J2 loOnm'n^Nnij). (28b) 

n=—cx) 

Eqs. (28a) and (28b) can be further simplified. Taking into account Eq. (26a) 
and the fact that l/{z + ie) = P(l/z) — m5{z) we can write 

where we took into account that D{uj) is an odd function. When summing over 
n in Eq. (28a), only the last term in Eq. (29) will contribute for every Since 
the contributions from the first term in D„ and D-n cancel for each n because 

D{uj) D{uj) _ 2nD{uj) 

Lo ~ nQ{j) to + nQ{j) iv"^ — n'^Q'^{j) 

is an odd function and integrates to zero. The final result for the coefficient V 
is 

oo 

'^{j)=i Yl K{j)\^nD{nn). (30a) 

n=— oo 

Performing similar steps with Nn, we get the following expression for A/": 

oo 

A^(i) = 2 E Mj)fnm{nn). (30b) 

n=— oo 

Particularizing to the case of an ohmic environment initially at zero temper- 
ature, we have N[(jj) = r]\ijj\, D{uj) = ir]LO, which lead to 

oo 

T>{j)=2nn{j)J2\xn{j)\^n\ (31a) 

n=0 

oo 

n=0 

4 Environment induced decay rate 

Our aim in this section is to compute the environment induced decay rate by 
solving the averaged Fokker-Planck equation, Eq. (27), with the appropriate 
boundary conditions. We begin by introducing a simplifying hypothesis. 



V 



d9 dx^ 

V 

27r de 



d9 dx dN 
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4.1 The hcirmonic approximation 

The non-dissipative dynamics described by (18) is approximately harmonic for 
energies much lower than the escape energy eg. In this case the frequency of 
that motion is simply Qq. For intermediate energies, the motion is qualitatively 
similar, but with a somewhat smaller frequency. It is only for energies extremely 
close to the escape energy eg that the motion is substantially different: the 
particle needs a very large amount of time to complete one period and hence 
the frequency tends to zero. 

In order to obtain an order of magnitude estimate of the escape rate, it 
is legitimate to approximate the classical motion of Eq. (18) by its harmonic 
approximation, since this approximation is qualitatively valid for all the values 
of the energy except for a very small region of energies extremely close to Eg-^ 

In the case of vacuum decay in quantum field theory studied by Calzetta 
et al. (2001, 2002) the harmonic approximation could not be introduced, since 
in that case there existed a frequency threshold with a value greater than 
in the dissipation, so that only Fourier modes with a frequency higher than 
the threshold contributed to it. Hence it was crucial to consider the fully non- 
linear dynamics of the system. On the other hand, in our case the noise and 
dissipation kernels do not exhibit such a threshold and it is possible to introduce 
the harmonic approximation, which amounts to neglect the Fourier modes at 
higher frequency in front of the lowest ones in the solution of Eq. (18). 

Neglecting the nonlinear term of the potential in Eq. (18b), the solution of 
the equations of motion with final conditions X{t) = x and P{t) = p can be 
written as 

P{t') = ^ (iMnox + p) e^^o(*'-*) +^ {-iM9.QX + p) e-^"o(*'-*) . (32b) 

The frequency of the motion is simply given by fi(j) = Hq, and thus the action 
variable j is simply given by j = e/Qq, being e = p^/(2M) + J7ga;^/2 the energy. 
In order to determine the action variable, we compare the solution for X{t') 
with Eq. (24), and write 

^C' = 4^1^ 

where we have identified the angle variable as 
V(MOox)2+j*2- 

^In fact, an exact solution of the equations of motion, together with a numerical analy- 
sis of the Fokker-Plank equation, reveal that the result that we will obtain here is not only 
qualitatively, but also quantitatively valid within the degree of approximation we are working. 
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The decomposition of the x variable in terms of 9 and j is given by 



1 / V .„ . 1 / 2i 



SO that 



=x_i(i) = ^y^J^, a;„ = 0, n/-l,l. (36) 
The coefficients T>{j) and N{j) can be therefore expressed as 

Vij)=2jj, Srij)=7j, (37) 

[we recall that 7 = r//(2M)], and the averaged Fokker-Planck equation may be 
written as 

dF ^ d fhflo dF ^ \ 

or, equivalently, working with energies, 

dF ^ d f dF \ 

where £0 = It can also be rewritten as a conservation equation, 

dF 

where 

^ = -2^ (^e^e^ + eF^ (41) 

is the probability flux. 

4.2 Escape rate: Normal mode analysis 

Assuming that the Fokker-Planck equation can be decomposed into a sum of 
normal modes, 

F{t,E)=Y,Cre-'-'fr{e), (42) 

r 

we get the following time-independent equation: 

Lfr + rfr = 0, L = 27-^ (eoe^ + e] . (43) 



de \ de 

The boundary conditions of the partial differential equation are the following. 
First, we have assumed that the particle is removed once it arrives at the sepa- 
ratrix, so that there will be no probability to find the particle at the separatrix: 
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fri^s) = 0- Secondly, we will also assume a vanishing flux of incoming particles 
at £ = 0, i.e., $(0) = 0, which is equivalent to demanding the finiteness of fr 
and its derivative at e = (/r(0), /^(O) < oo), as can be seen from Eq. (41) . 

The normal-mode analysis can be formulated as a standard Sturm-Liouville 
problem. The operator L, which may be written as 

L = 27 e-^/^° (e'/'<^ e-^^ + 27, (44) 

is self-adjoint with the aforementioned boundary conditions and the scalar prod- 
uct defined by 

if, 9) = (27)-^ r d£ e-/^o ne)gie). (45) 
Jo 

Thus, the theory of differential equations guarantees that the eigenf unctions 

/,.(£) constitute a complete orthogonal set, and that the eigenvalues r are real 
(Courant and Hilbert, 1953). Furthermore, the operator L is negative definite, 
which can be seen as follows: 

{f,Lf) = £0 p d£ e^/^o /*(£)^ e e-^/^o A (^e^Ao 



'd£ £6-^/^0 



d • 



d£ 



(e^/^o/(£)) 



< 0, / / 0, 



where we have integrated by parts in the last equality. This implies that the 
eigenvalues r are always positive, as expected. 

Performing the change of variables fr{s) = 27 e"^' ffiu), where y = e/eo and 
f = r/(27), the differential equation Lfr = rfr adopts the form of the Laguerre 
differential equation, 

yffiy) + (1 - y)/^(y) + f/^y) = 0, (46) 
whose unique regular solution is given by 

My) = NLr-{y), (47) 

where Lr{y) are Laguerre functions (which reduce to the Laguerre polynomials 
in the case of non- negative integer f; see Gradsteyn and Ryzhik, 1980), and N 
is a normalization constant. 

The solution we have found verifies the first of the boundary conditions, the 
regularity at the origin. Now we impose the second of the boundary conditions, 
namely ff{ys) = (with y^ = £s/eo)- This boundary condition will imply a 
discretization on the possible values of the escape rate r. Finally, knowing 
the initial state /i(£) = F(0, £), we will be able to reconstruct the solution by 
computing the coefficients c^: 



(48) 
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Unfortunately, it is not possible to determine analytically the possible values 
of f from the equation Lf (ys) = 0- However, if the potential barrier were in- 
finitely high, the second boundary condition would read limj^^oo Lf{y) = 0, 
implying that Lf{y) should have, at most, a polynomial behavior at infinity, and 
the eigenvalues would be f = 0, 1, 2 . . ., so that eigenmodes with f ^ would 
decay in a time given by 7"^ or shorter. However, the potential barrier corre- 
sponds to a some large but finite energy, and hence the real eigenvalues differ 
from those computed in the infinite barrier case by a small quantity, at least for 
those eigenvalues corresponding to eigenstates with characteristic energies much 
lower than the potential barrier. Thus, wc can compute them perturbatively. 

Furthermore, the relevant contribution to the decay rate will be given by the 
mode with the lowest eigenvalue, which will fulfill the condition f <^ 1, since the 
remaining modes will decay in a time of order 7"^ at most. Thus, we proceed 
to compute perturbatively the lowest order mode by expanding the Laguerre 
function around f = 0: 

Lr{y) = 1 + f[lny + 7^ - Ei(y)] + 0(f2), (49) 

where 7b = 0.577216 ... is Euler's constant, and Ei(?/) is the exponential integral 
defined as Ei(?/) = PV Jl^{^^/u) du. The expansion in Eq. (49) can be found by 
solving perturbatively Eq. (46) up to order f , and imposing the correct boundary 
condition at y = 0. 

Therefore, imposing the boundary condition Lf{ys) = is equivalent to 
demanding 

r = 2jf^———- — , , . (50) 

Ei(es/eo) - ln(es/eo) - 7e ^ ' 

Since Eg is much larger than eo, the exponential integral can be approximated by 
Ei(es/£o) ~ (eo/^s) 6^=/^°, and the other two terms in the denominator become 
negligible in front of this one. Hence, we may approximate the lowest order 
solution by 

■ exp . (51) 



^0 V ^0 , 

Equation (51), which gives the probability per unit time for a particle to jump 
the barrier, is our final result for the escape rate. 

In the case of vacuum decay in quantum field theory studied by Calzetta et al. 
(2001, 2002) the time-independent Fokker-Planck equation had a continuous 
spectrum, and therefore in that case it was not possible to perform either an 
eigenvalue expansion, like the one we have performed in this subsection, or follow 
Kramers's method, as will be done in the next subsection. 
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4.3 Escape rate: Kramers's method 

Our Fokker-Planck equation (39) is completely analogous to that found by 
Kramers (1940) in the classical underdamped case, once we replace kT, the tem- 
perature times Boltzmann's constant, by the zero point energy of the harmonic 
oscillator eq = hQ.o/2? Therefore, as an alternative to the previous subsection, 
we can apply the same method as Kramers in order to compute the decay rate. 

Instead of imposing the correct boundary conditions, we look for the solu- 
tions of the Fokker-Planck equation with constant flux $ = $o- Since the flux 
can be rewritten as 

$ = -27 e-^/^o eeo-^ (^e'/'° , (52) 
these solutions will be given by 



F(e) = e-=/=o / de' 

^ ' 27£o Je 



^0 

27^0 



(53) 



[Ei(es/£o) - Ei(£/eo)] • 



Again, we have imposed that -P(es) = because we assume that when a particle 
arrives at the scparatrix it never reenters the potential well region. Notice that 
the solution we have found has a logarithmic singular behavior at small energies, 
which takes into account the injection of a probability flux through the point 
£ = necessary for the maintenance of the constant flux. However, since we 
expect the flux of probability to be very small, this contribution will be not very 
significative and will affect only the region of energies e ^ Eq, which are much 
lower than the scape energy Es- 

We can compute the flux $0 by imposing the correct normalization of the 
averaged Wigner function F{e), namely jQ^d£F{E) = 1: 

/ d£ / de' —— = 1. (54) 

The main contribution to the integral over sf is due to those values for e' which 
differ from by a quantity of order eq. We may replace e' by its value at 
the separatrix, e^- Making this approximation we can perform analytically the 
integrals in Eq. (54). Retaining only those terms which have an exponential 
factor e^«/^o, we obtain the final result for the flux of particles: 

$ = $o~Mexpf-^V (55) 
£0 \ £oJ 



^Notice that this is only true under the harmonic approximation. Had not we neglected the 
cubic term in the equations of motion, our Fokker-Planck equation would not be equivalent to 
the one found by Kramers in the classical underdamped case. 
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Although we have assumed that the flux is constant in order to solve the differ- 
ential equation, this is not actually the case, and the flux at the separatrix <l>s is 
proportional to the total probability of finding the particle in the potential well 
region: ^s{t) = P(t)^o. Integrating the conservation equation (40) with respect 
to the energy, we easily find that the total probability decay follows the law 

which may be integrated to give P{t) = e^*"*. Hence the probability flux $0 
is to be identified with the escape rate r of last subsection. We see that both 
methods agree. 



5 Discussion 

The study of quantum tunneling based on a real time formulation seems crucial 
in order to address highly non-equilibrium situations in which no adiabaticity 
assumptions can be made. As mentioned in the introduction, a first step in this 
direction was made in Calzetta et al. (2001, 2002), where the effect on vacuum 
decay in quantum field theory due to the backreaction of the short-wavelength 
modes was analyzed by regarding the long-wavelength modes responsible for 
tunneling as an open quantum system. In that case it was found that such a 
backreaction seemed to yield an enhancement of the decay rate. There are, how- 
ever, a couple of aspects which deserve, in our opinion, a more careful analysis. 
The first one is the need for a suitable identification of the tunneling degrees 
of freedom and the corresponding implementation of a system-environment sep- 
aration which leads to the reformulation of the problem in terms of an open 
quantum system. Secondly, when solving the master equation that governs the 
time evolution of the reduced Wigner function for the tunneling degree of free- 
dom in order to obtain the vacuum decay rate, the attention was focused on the 
backreaction of the environment (the short-wavelength modes) and the higher 
derivative terms which would be uniquely responsible for tunneling if the system 
were isolated were neglected. 

In this paper we have considered a fairly simple quantum mechanical open 
system, in which the system-environment separation is given before-hand and 
the coupling constant governing the interaction between the system and the 
environment can be adjusted at will, rather than being sclf-consistently deter- 
mined, as happened to be the case for vacuum decay in field theory. Therefore, 
it is interesting to check whether a similar enhancement of the tunneling rate 
is obtained in that simpler model and see if such an effect is generic. Further- 
more, the result can be considered more robust than that in the field theory case 
since the number of important assumptions made is rather small. It is, thus. 
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worthwhile to elaborate on this point and recall the different approximations 
employed throughout the paper in order to obtain the decay rate. 

First of all, the master equation that we are considering in this paper can be 
obtained by keeping terms of order h'^ and h'y at most, and restricting to values 
of 7 which are small enough (in particular we should take (7/^0)5' < h, where 
S is the characteristic action of the problem) , so that the terms of order 7^ or 
higher can be neglected^ (Roura, 2001). On the other hand, the master equation 
obtained in this way contains a term with third order derivatives of the reduced 
Wigner function with respect to the momentum (the term of order h^). This sort 
of terms, which arc absent in any diffusion equation with a finite Kramers-Moyal 
expansion associated to a classical stochastic process, are intimately related to 
genuinely quantum effects due to the nonlinearities of the potential and imply 
that even a reduced Wigner function which is initially positive everywhere will 
acquire negative values when it evolves in time. Moreover, this higher derivative 
term would be uniquely responsible for tunneling if the system where isolated. 
Despite its remarkable features, this term has been neglected in the present 
article, since otherwise we were unable to obtain analytical results for the decay 
rate. This approximation, whose justification will be further discussed below, 
constitutes the most drastic assumption made throughout the paper. 

Having neglected the higher derivative terms, the master equation is equiv- 
alent to the Fokker-Planck equation associated to a Langevin equation with a 
Gaussian stochastic source characterized by a nonlocal correlation function (the 
nonlocal noise kernel). From this point on, most of the approximations em- 
ployed in order to compute the decay rate are more or less standard (Hanggi 
et al, 1990). First, we change to action-angle variables and make use of an adia- 
batic approximation to eliminate the fast variable (the angle) . This is consistent 
provided that 7 ^ f^O) which is in agreement with the previous assumption of 
small enough values for 7. Next, a harmonic approximation is introduced for 
the solutions to the equations of motion for the isolated system which appear in 
both the master equation and the Fokker-Planck equation. This approximation 
helps to obtain a rather simple result for the decay rate and can be justified both 
qualitatively and quantitatively, in contrast to the field theory case analyzed by 
Calzetta et al. (2001, 2002), where the existence of a threshold for the dissi- 
pation and noise kernels would preclude such an approximation. Finally, it is 
assumed that the characteristic decay time is much smaller than the relaxation 
and dynamical timescales: r ^ 7 <C i^Oj where r is the decay rate (the inverse 

■^One might bo concerned that the truncation of higher orders in 7 when considering times 
much larger than the characteristic relaxation timescale (i.e., yt 3> 1), as required in order to 
compute the decay rate, might no longer be valid due to the existence of secular terms among 
the terms of higher order in 7 that have been neglected. Although arguments that justify such 
a truncation when computing the decay rate can be given, this point might deserve a more 
careful analysis. 
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of the decay time). Although 7 was required to be small enough, the previous 
inequality can be fulfilled provided that the potential barrier is sufficiently high, 
i.e., Ss is large enough. 

After reviewing the main approximations employed, let us compare the result 
obtained for the decay rate to that of Caldeira and Leggett (1981, 1983b) as 
well as to the tunneling rate when the coupling to the environment is absent. 
Whereas Caldeira and Leggett found that the interaction with the environment 
tends to suppress tunneling, we are not going to recover that result since, having 
neglected the higher derivative terms responsible for tunneling in an isolated 
system, the dissipation and diffusion terms that appear in the master equation 
due to the backrcaction of the environment will always lead to a positive (or, at 
most, vanishing) probability of escaping from the metastable well. Nevertheless, 
if the decay rate obtained were much smaller or much larger than the tunneling 
rate for the isolated system so that the timescales governing both processes are 
very different, one could expect that the contribution to the decay rate from the 
process with a shortest characteristic timescale would be dominant. 

The tunneling rate for an isolated system initially trapped in the metastable 
minimum of the potential considered in this paper is (Caldeira and Leggett, 
1983b) 

rt ~ rio — exp — , (56) 



,£0/ V 5 £0, 

where we recall that eo is the zero-point energy of a harmonic oscillator of fre- 
quency ^0- On the other hand, Caldeira and Leggett (1981, 1983a) obtained the 

modification of the tunneling rate due to the interaction with the environment, 
which in the case of small dissipation is given by 



1/2 

rt ~ r^o ( — ) exp 



£s /^18^54C(3) 7 



£0 \ 5 rii 



(57) 



Therefore, since the interaction with the environment simply adds a negative 

contribution to the exponent, it always tends to suppress the tunneling rate. 
Finally, the decay rate due to the activation-like effect obtained in the previous 
section corresponds to 

r ^ 7— exp ( — ^ J , (58) 



Co V ^0 , 

which is valid for small 7 and large Eg ■ When 7 is very small the decay rates from 
Eqs. (56) and (57) become very close. Furthermore, although a small 7 implies 
a small contribution to the decay rate associated to the activation-like process, 
it can always be made arbitrarily larger than the tunneling rate by taking Sg 
large enough. 

The fact that the activation-like decay rate can be made arbitrarily large 
as compared to the tunneling rate for the isolated system seems to suggest, as 
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mentioned above, that even if the two phenomena were considered simultane- 
ously, the former would be expected to dominate. If that were the case, it would 
imply that the usual instanton methods, when applied to a system interacting 
with an environment, downplay the role of the backreaction of the environment 
on the system dynamics. This can be qualitatively interpreted in the following 
way: while, roughly speaking, the tunneling effect for an isolated system can 
be regarded as a consequence of the energy fluctuations implied by Heisenberg's 
uncertainty principle, the interaction with an environment would induce fluctu- 
ations on the system due to the quantum fluctuations of the environment itself 
(Nagaev and Biittiker, 2002), which would enhance the tunneling rate. 

Nevertheless, in order to reach a definite conclusion it would be desirable to 
deal with the two contributions simultaneously and make sure that the higher 
derivative terms do not disrupt the effect of the backreaction terms, even when 
the timescale for the contribution to the decay rate from the former terms is 
much longer than the timescale associated to the activation-like process. Un- 
fortunately, dealing with the higher derivative terms is not an easy task, and it 
seems hard to provide a real time description of the tunneling process in terms of 
the Wigner function even for an isolated system; see, however, Riskcn and Vogel 
(1988) for a first step in that direction. One may try to gain partial information 
on the relation between the two processes by considering different potentials 
with wider or narrower barriers, since one would naively expect that tunnel- 
ing is suppressed for wide barriers while the activation-like contribution does 
not depend so much on the width of the barrier, as long as the height remains 
the same. A couple of comments concerning the freedom to modify the system 
potential are, nevertheless, in order. First, the potential must be analytic in 
order the derive the master equation for the reduced Wigner function. Second, 
when solving the formally equivalent problem of a classical Brownian particle 
escaping from the potential well, we assumed that once the particle reaches the 
maximum of the barrier it escapes and never comes back, but if a very wide 
barrier is considered, the probability that the particle comes back due to the 
fluctuations may become no longer negligible. 

We close this section insisting on the importance of finding a satisfactory 
method to deal with the higher derivative terms, which would be very helpful in 
order to elucidate whether the enhancement of the decay rate obtained in this 
paper and entirely due to the back reaction from the environment fluctuations 
would still persist when the terms responsible for tunneling in isolated systems 
are also taken into account. Such a method would have an interest in its own 
right even if the results of Caldeira and Leggett (1981, 1983b) were finally re- 
covered when properly taking into account the higher derivative terms, since it 
would constitute a key step in formulating a real time description of tunneling. 
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